###1.simualtion from different distributions:

##normal distribution:
rnorm(10,mean=0,sd=2)

##chi-square distribution:
rchisq(100, df=2) 

##exponential distribution:
rexp(100, 2)

##t distribution:
rt(10,df=2) 

##Uniform distribution:
runif(10, min=1, max=2)

##Gamma distribution:
rgamma(20, shape=1, rate=2)

##Bernoulli distribution:
library(Rlab)
rbern(10, 0.5)

##Binomial distribution:
rbinom(10, 2, prob=0.5)  

##Poisson distribution:
rpois(20,1)  

##Multivariate normal distribution:
library(MASS)
mvrnorm(n=10, mu=c(1,2), Sigma=matrix(c(1,0,0,1),2,2))

##Simple random samples:
sample(1:10, 2, replace = FALSE)
sample(1:10, 5, replace = TRUE)

##Example 1: One sample problem:

##mod=1 for normal distribution, mod=2 for chi-square distribution
##B: Monte Carlo sample size
##n: sample size for each Monte Carlo sample

fsim1<-function(mod,B,n){
Etheta<-NULL
EV<-NULL
EsI<-NULL
EL<-NULL
iter<-0
repeat{
iter<-iter+1
if(mod==1){
           Y<-rnorm(n,mean=1,sd=1)
           theta0<-1
          }
else if(mod==2){
                Y<-rchisq(n,df=2)
                theta0<-2
               }
etheta<-mean(Y)
ev<-1/n*var(Y)
LB<-etheta-1.96*sqrt(ev)
UB<-etheta+1.96*sqrt(ev)
sI<-as.numeric(theta0>LB & theta0<UB)
eL<-UB-LB
Etheta<-c(Etheta,etheta)
EV<-c(EV,ev)
EsI<-c(EsI,sI)
EL<-c(EL,eL)
if(iter==B){break}
} 

BIAS<-mean(Etheta)-theta0
SE<-sqrt(var(Etheta))
MSE<-BIAS^{2}+SE^{2}
RB<-BIAS/theta0
RSE<-SE/theta0
RRMSE<-sqrt(MSE)/theta0
VRB<-(mean(EV)-SE^{2})/SE^{2}
CR<-mean(EsI)
AL<-mean(EL)

res<-c(BIAS,SE,MSE,RB,RSE,RRMSE,VRB,CR,AL)
return(res)
}

fsim1(1,2000,200)
fsim1(2,2000,200)

##Example 2: Two sample problem:

##B: Monte Carlo sample size
##n1: sample size of sample one for each Monte Carlo sample
##n2: sample size of sample two for each Monte Carlo sample

fsim2<-function(B,n1,n2){
Etheta<-NULL
EV<-NULL
EsI<-NULL
EL<-NULL
iter<-0
repeat{
iter<-iter+1
Y1<-rnorm(n1,mean=2,sd=1)
Y2<-rnorm(n2,mean=1,sd=1)

theta0<-2-1
          
etheta<-mean(Y1)-mean(Y2)
ev<-(1/n1+1/n2)*((n1-1)*var(Y1)+(n2-1)*var(Y2))/(n1+n2-2)
LB<-etheta-1.96*sqrt(ev)
UB<-etheta+1.96*sqrt(ev)
sI<-as.numeric(theta0>LB & theta0<UB)
eL<-UB-LB
Etheta<-c(Etheta,etheta)
EV<-c(EV,ev)
EsI<-c(EsI,sI)
EL<-c(EL,eL)
if(iter==B){break}
} 

BIAS<-mean(Etheta)-theta0
SE<-sqrt(var(Etheta))
MSE<-BIAS^{2}+SE^{2}
RB<-BIAS/theta0
RSE<-SE/theta0
RRMSE<-sqrt(MSE)/theta0
VRB<-(mean(EV)-SE^{2})/SE^{2}
CR<-mean(EsI)
AL<-mean(EL)

res<-c(BIAS,SE,MSE,RB,RSE,RRMSE,VRB,CR,AL)
return(res)
}

fsim2(2000,200,200)

##Example 3: Linear regression:

##B: Monte Carlo sample size
##n: sample size for each Monte Carlo sample

fsim3<-function(B,n){
Etheta<-NULL
EV<-NULL
EsI<-NULL
EL<-NULL
iter<-0
repeat{
iter<-iter+1
X<-rnorm(n,mean=0,sd=1)
epsilon<-rnorm(n,mean=0,sd=1)
Y<-1+X+epsilon
theta0<-c(1,1)
fit<-lm(Y~X)
etheta<-summary(fit)$coefficients[,1]
ev<-(summary(fit)$coefficients[,2])^{2}

LB<-etheta-1.96*sqrt(ev)
UB<-etheta+1.96*sqrt(ev)
sI<-as.numeric(theta0>LB & theta0<UB)
eL<-UB-LB
Etheta<-cbind(Etheta,etheta)
EV<-cbind(EV,ev)
EsI<-cbind(EsI,sI)
EL<-cbind(EL,eL)
if(iter==B){break}
} 

BIAS<-apply(Etheta,1,mean)-theta0
SE<-sqrt(apply(Etheta,1,var))
MSE<-BIAS^{2}+SE^{2}
RB<-BIAS/theta0
RSE<-SE/theta0
RRMSE<-sqrt(MSE)/theta0
VRB<-(apply(EV,1,mean)-SE^{2})/SE^{2}
CR<-apply(EsI,1,mean)
AL<-apply(EL,1,mean)

res<-cbind(BIAS,SE,MSE,RB,RSE,RRMSE,VRB,CR,AL)
return(res)
}

fsim3(2000,200)

##Example 4: Logistic regression:

##B: Monte Carlo sample size
##n: sample size for each Monte Carlo sample
##M: Bootstrap sample size

fsim4<-function(B,n,M){
Etheta<-NULL
EV1<-NULL
EsI1<-NULL
EL1<-NULL

EV2<-NULL
EsI2<-NULL
EL2<-NULL

iter<-0
repeat{
iter<-iter+1
X<-rnorm(n,mean=0,sd=1)
p<-exp(1+X)/(1+exp(1+X))
Y<-rbern(n,p)
theta0<-c(1,1)
fit<-glm(Y~X,family=binomial(link = "logit"))
etheta<-summary(fit)$coefficients[,1]
ev1<-(summary(fit)$coefficients[,2])^{2}

##Bootstrap variance estimator:
DAT<-cbind(Y,X)
b<-0
Btheta<-NULL
repeat{
b<-b+1
b_ID<-sample(1:n,n,replace=TRUE)
bDAT<-DAT[b_ID,]
bfit<-glm(bDAT[,1]~bDAT[,2],family=binomial(link = "logit"))
btheta<-summary(bfit)$coefficients[,1]
Btheta<-cbind(Btheta,btheta)
if(b==M){break}
}
ev2<-apply(Btheta,1,var)

LB1<-etheta-1.96*sqrt(ev1)
UB1<-etheta+1.96*sqrt(ev1)
sI1<-as.numeric(theta0>LB1 & theta0<UB1)
eL1<-UB1-LB1

LB2<-etheta-1.96*sqrt(ev2)
UB2<-etheta+1.96*sqrt(ev2)
sI2<-as.numeric(theta0>LB2 & theta0<UB2)
eL2<-UB2-LB2

Etheta<-cbind(Etheta,etheta)
EV1<-cbind(EV1,ev1)
EsI1<-cbind(EsI1,sI1)
EL1<-cbind(EL1,eL1)

EV2<-cbind(EV2,ev2)
EsI2<-cbind(EsI2,sI2)
EL2<-cbind(EL2,eL2)

if(iter==B){break}
} 

BIAS<-apply(Etheta,1,mean)-theta0
SE<-sqrt(apply(Etheta,1,var))
MSE<-BIAS^{2}+SE^{2}
RB<-BIAS/theta0
RSE<-SE/theta0
RRMSE<-sqrt(MSE)/theta0
VRB1<-(apply(EV1,1,mean)-SE^{2})/SE^{2}
CR1<-apply(EsI1,1,mean)
AL1<-apply(EL1,1,mean)

VRB2<-(apply(EV2,1,mean)-SE^{2})/SE^{2}
CR2<-apply(EsI2,1,mean)
AL2<-apply(EL2,1,mean)

res<-cbind(BIAS,SE,MSE,RB,RSE,RRMSE,VRB1,CR1,AL1,VRB2,CR2,AL2)
return(res)
}

fsim4(1000,100,100)

> fsim4(1000,100,100)
                  BIAS        SE        MSE         RB       RSE     RRMSE
(Intercept) 0.02586266 0.2597998 0.06816481 0.02586266 0.2597998 0.2610839
X           0.05523762 0.2972840 0.09142899 0.05523762 0.2972840 0.3023723
                    VRB1   CR1      AL1      VRB2   CR2      AL2
(Intercept)  0.003757506 0.962 1.015219 0.1487553 0.967 1.077950
X           -0.005273640 0.960 1.150720 0.1843197 0.967 1.241654